package scu.maqiang.fes;

import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES1L2PE extends FES {
	public FES1L2PE(Mesh1L2 mesh1D) {
		super(mesh1D, 2);
		this.eleStiff = FES1L2PE::elementPEStiff;
		this.eleMass = FES1L2PE::elementPEMass;
		this.eleSource = FES1L2PE::elementPESource;
	}
	
	/**
	 * 
	 * @param coefFunc
	 * @param param
	 * @param tp
	 * @param A
	 */
	public void assemblePEStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES1L2PE::elementPEStiff, tp, A);
	}
	
	public void assemblePEStiff(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES1L2PE::elementPEStiff, tp, A);
	}

	public void assemblePEStiff(double[][] vecCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(vecCoef, FES1L2PE::elementPEStiff, tp, A);
	}
	
	public static void elementPEStiff(int[] ele, double[][] xx, double[][] coef, BVPType tp, double[][] eM) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[0] - x[1]);
		double mx = 0.5 * (x[0] + x[1]);
		//double[] dx = { -1.0 / L, 1.0 / L};
		//double[][] mat = new double[][] {{mG, mEp1}, {mEp1, -mE1}};
		if (tp == BVPType.COMMON) {
			double mG = 0.5 * (coef[0][0] + coef[0][1]);
			double mEp = 0.5 * (coef[1][0] + coef[1][1]);
			double mE1 = 0.5 * (coef[2][0] + coef[2][1]);
			eM[0][0] = eM[2][2] = mG / L;
			eM[0][2] = eM[2][0] = -mG / L;
//			eM[0][1] = eM[1][0] = eM[2][3] = eM[3][2] = mE1 / L;
//			eM[0][3] = eM[3][0] = eM[1][2] = eM[2][1] = -mE1 / L;
//			eM[1][1] = eM[3][3] = -mEp / L;
//			eM[1][3] = eM[3][1] = mEp / L;
			eM[0][1] = eM[1][2] = eM[2][3] = eM[3][0] = mE1 / L;
			eM[1][0] = eM[2][1] = eM[3][2] = eM[0][3] = -mE1 / L;
			eM[1][1] = eM[3][3] = mEp / L;
			eM[1][3] = eM[3][1] = -mEp / L;
			//System.out.println(MVO.toString(eM));
			
		} else if (tp == BVPType.AXISYMMETRIC){
			double mG1 = 0.5 * (coef[0][0] + coef[0][1]);
			double mG2 = 0.5 * (coef[1][0] + coef[1][1]);
			double lambda = 0.5 * (coef[2][0] + coef[2][1]);
			double mEp = 0.5 * (coef[3][0] + coef[3][1]);
			double mE1 = 0.5 * (coef[4][0] + coef[4][1]);
			double mE2 = 0.5 * (coef[5][0] + coef[5][1]);
			double[][] D = {{mG1, lambda, mE1}, {lambda, mG2, mE2}, {-mE1, -mE2, mEp}};
			double[][] BT = new double[4][3];
			double[][] B = new double[3][4];
			BT[0][0] = BT[1][2] = -1.0 / L;
			BT[0][1] = BT[2][1] = 1.0/(x[0] + x[1]);
			BT[2][0] = BT[3][2] = 1.0 / L;
			
			B[0][0] = B[2][1] = -mx;
			B[1][0] = B[1][2] = 0.5 * L;
			B[0][2] = B[2][3] = mx;
			MVO.matmul(BT, MVO.matmul(D, B), eM);
		} else if (tp == BVPType.CONSTITUTE_MATRIX_COMMON) {
			
		} else if (tp == BVPType.CONSTITUTE_MATRIX_AXIS) {
			
		}
		else {	throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assemblePESource(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES1L2PE::elementPESource, tp, RHS);
	}

	public void assemblePESource(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES1L2PE::elementPESource, tp, RHS);
	}

	public void assemblePESource(double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(vecCoef, FES1L2PE::elementPESource, tp, RHS);
	}

	public static void elementPESource(int[] ele, double[][] xx, double[][] vf, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] f = vf[0];
		double[] rf = vf[1];
		if (tp == BVPType.COMMON) {
			eleF[0] = (2 * f[0] + f[1]) * L / 6;
			eleF[1] = (2 * rf[0] + rf[1]) * L / 6;
			eleF[2] = (f[0] + 2 * f[1]) * L / 6;
			eleF[3] = (rf[0] + 2 * rf[1]) * L / 6;
		} else if (tp == BVPType.AXISYMMETRIC) {
			double sxsf = MVO.sum(x) * MVO.sum(f);
			double sxsrf = MVO.sum(x) * MVO.sum(rf);
			eleF[0] = (2 * x[0] * f[0] + sxsf) * L / 12.0;
			eleF[1] = (2 * x[0] * rf[0] + sxsrf) * L / 12.0;
			eleF[2] = (2 * x[1] * f[1] + sxsf) * L / 12.0;
			eleF[3] = (2 * x[1] * rf[1] + sxsrf) * L / 12.0;
		} else if (tp == BVPType.SPHERICAL) {
			
		} else {
			throw new IllegalArgumentException("Wrong Type!");
		}
	}
	
	public void assemblePEMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES1L2PE::elementPEMass, tp, A);
	}

	public void assemblePEMass(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES1L2PE::elementPEMass, tp, A);
	}

	public void assemblePEMass(double[][] vecCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(vecCoef, FES1L2PE::elementPEMass, tp, A);
	}
	
	public static void elementPEMass(int[] ele, double[][] coord, double[][] g, BVPType tp, double[][] k) {
	    double[] coefM = g[0];
	    double[] x = MVO.col(coord, 0);
	    double L = Math.abs(x[0] - x[1]);
	    MVO.fill(k, 0.0);
	    switch(tp) {
	    case CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] + coefM[1]) * L;
	    	k[0][2] = k[2][0] = (coefM[0] + coefM[1]) * L;
	        k[2][2] = (coefM[0] + 3 * coefM[1]) * L;
	        break;
	    case LUMP_MASS:
	    	L /= 6.0;
	    	k[0][0] = (2 * coefM[0] + coefM[1]) * L;
	    	k[2][2] = (coefM[0] + 2 * coefM[1]) * L;
	    	break;
	    case AXIS_CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] * x[0] + coefM[1] * x[1]) * L;
	    	k[0][2] = k[2][0] = (coefM[0] * x[0] + coefM[1] * x[1]) * L;
	        k[2][2] = (coefM[0] * x[0] + 3 * coefM[1] * x[1]) * L;
	        break;
	    case AXIS_LUMP_MASS:
	        L /= 6.0;
	        k[0][0] = (2 * coefM[0] * x[0] + coefM[1] * x[1]) * L;
	        k[2][2] = (coefM[0] * x[0] + 2 * coefM[1] * x[1]) * L;
	        break;
	    case SPHE_CONSIST_MASS:
	    	L /= 12.0;
	    	k[0][0] = (3 * coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	    	k[0][2] = k[2][0] = (coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	        k[2][2] = (coefM[0] * x[0] * x[0] + 3 * coefM[1] * x[1] * x[1]) * L;
	        break;
	    case SPHE_LUMP_MASS:
	        L /= 6.0;
	        k[0][0] = (2 * coefM[0] * x[0] * x[0] + coefM[1] * x[1] * x[1]) * L;
	        k[2][2] = (coefM[0] * x[0] * x[0] + 2 * coefM[1] * x[1] * x[1]) * L;
	        break;
	    default:
	    	throw new IllegalArgumentException("Wrong BVPType!");
	    }	   
	}
	
	
	public void assemblePressure(ScalarFunc[] coeffFunc, double[] param, BVPType tp, double[] RHS,
			int... bdLabel) {
		assembleBoundaryVector(coeffFunc, param, FES1L2PE::elementPressure, tp, RHS, bdLabel);
	}
	
	public void assemblePressure(double[] constCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(constCoef, FES1L2PE::elementPressure, tp, RHS, bdLabel);
	}

	public void assemblePressure(double[][] vecCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(vecCoef, FES1L2PE::elementPressure, tp, RHS, bdLabel);
	}
	
	public static void elementPressure(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] k) {
		double x = coord[0][0];
		switch (tp) {
		case COMMON:
			k[0] = g[0][0];
			break;
		case AXIS_PLANESTRAIN:
		case AXIS_PLANESTRESS:
		case AXISYMMETRIC:
			k[0] = g[0][0] * x;
			break;
		case SPHERICAL:
			k[0] = g[0][0] * x * x;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void computeGradient(double[] u, double[][] gradient) {
	    double[] count = new double[ndof];
	    double[][] coord = new double[2][1];
	    double[] eleU = new double[2];
	    double[] elePhi = new double[2];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	int[] dof = getDof(ele);
	    	eleU[0] = u[dof[0]];
	    	eleU[1] = u[dof[2]];
	    	elePhi[0] = u[dof[1]];
	    	elePhi[1] = u[dof[3]];
	    	
	    	gradient[0][dof[0]] += eleU[1] - eleU[0];
	    	gradient[0][dof[1]] += elePhi[1] - elePhi[0];
	    	gradient[0][dof[2]] += eleU[1] - eleU[0];
	    	gradient[0][dof[3]] += elePhi[1] - elePhi[0];
	    	
	    	count[ele[0]] += x[1] - x[0];
	    	count[ele[1]] += x[1] - x[0];
	    }
	    for(int i = 0; i < ndof; i++) {
	    	gradient[0][i] /= count[i];
	    }
	    
	    double[] xx = { mesh.getNode(0)[0], mesh.getNode(1)[0], mesh.getNode(2)[0]};
	    //�ֱ����߽��ϵĵ���dU(1)��dU(N)
	    double dx1 = xx[1] - xx[0]; 
	    double dx2 = xx[2] - xx[1];
	    gradient[0][0] = ((u[2] - u[0]) * (dx1 + dx2) * (dx1 + dx2) 
	    		        - (u[4] - u[0]) * dx1 * dx1) / (dx1 * dx2 * (dx1 + dx2));
	    gradient[0][1] = ((u[3] - u[1]) * (dx1 + dx2) * (dx1 + dx2) 
		        - (u[5] - u[1]) * dx1 * dx1) / (dx1 * dx2 * (dx1 + dx2));
	    
	    int N = mesh.getNv();
	    xx[0] = mesh.getNode(N - 3)[0];
	    xx[1] = mesh.getNode(N - 2)[0];
	    xx[2] = mesh.getNode(N - 1)[0];
	    //dx1 = x(N) - x(N-1); dx2 = x(N - 1) - x(N - 2)
	    dx1 = xx[2] - xx[1]; dx2 = xx[1] - xx[0];
	    gradient[0][2 * N - 2] = -((u[2 * N - 4] - u[2 * N - 2]) * (dx1 + dx2)*(dx1 + dx2)  
	    		             - (u[2 * N - 6] - u[2 * N - 2]) * dx1*dx1) / (dx1 * dx2 * (dx1 + dx2));
	    gradient[0][2 * N - 1] = -((u[2 * N - 3] - u[2 * N - 1]) * (dx1 + dx2)*(dx1 + dx2)  
	             - (u[2 * N - 5] - u[2 * N - 1]) * dx1*dx1) / (dx1 * dx2 * (dx1 + dx2));
	}
	
	public void computeStressAndElectricDisp(MatrixFunc coefFunc, double[][] uPhi, double[][] ESED) {
	    double[] count = new double[ndof];
	    double[][] coord = new double[2][1];
	    double[] eleU = new double[2];
	    double[] elePhi = new double[2];
	    MVO.fill(ESED, 0.0);
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	int label = mesh.getElementLabel(i);
	    	double[] ip = mesh.getinteriorPointInElement(coord);
	    	double[][] D = coefFunc.action(ip, label, null);
	    	double[] x = MVO.col(coord, 0);
	    	double len = Math.abs(x[0] - x[1]);
	    	eleU[0] = uPhi[0][ele[0]];
	    	eleU[1] = uPhi[0][ele[1]];
	    	elePhi[0] = uPhi[1][ele[0]];
	    	elePhi[1] = uPhi[1][ele[1]];
	    	
 	    	ESED[0][ele[0]] += D[0][0] * (eleU[1] - eleU[0])/len + D[0][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[0][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[0][ele[1]] += D[0][0] * (eleU[1] - eleU[0])/len + D[0][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[0][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[1][ele[0]] += D[1][0] * (eleU[1] - eleU[0])/len + D[1][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[1][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[1][ele[1]] += D[1][0] * (eleU[1] - eleU[0])/len + D[1][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[1][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[2][ele[0]] += -D[2][0] * (eleU[1] - eleU[0])/len - D[2][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) - D[2][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[2][ele[1]] += -D[2][0] * (eleU[1] - eleU[0])/len - D[2][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) - D[2][2] * (elePhi[1] - elePhi[0]) / len;
	    	
	    	count[ele[0]] += 1.0;
	    	count[ele[1]] += 1.0;
	    }
	    
	    for(int i = 0, n = mesh.getNv(); i < n; i++) {
	    	ESED[0][i] /= count[i];
	    	ESED[1][i] /= count[i];
	    	ESED[2][i] /= count[i];
	    }
	}
	
	public void computeStressAndElectricDisp(double[][] D, double[][] uPhi, double[][] ESED) {
	    double[] count = new double[ndof];
	    double[][] coord = new double[2][1];
	    double[] eleU = new double[2];
	    double[] elePhi = new double[2];
	    MVO.fill(ESED, 0.0);
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	double len = Math.abs(x[0] - x[1]);
	    	eleU[0] = uPhi[0][ele[0]];
	    	eleU[1] = uPhi[0][ele[1]];
	    	elePhi[0] = uPhi[1][ele[0]];
	    	elePhi[1] = uPhi[1][ele[1]];
	    	
 	    	ESED[0][ele[0]] += D[0][0] * (eleU[1] - eleU[0])/len + D[0][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[0][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[0][ele[1]] += D[0][0] * (eleU[1] - eleU[0])/len + D[0][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[0][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[1][ele[0]] += D[1][0] * (eleU[1] - eleU[0])/len + D[1][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[1][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[1][ele[1]] += D[1][0] * (eleU[1] - eleU[0])/len + D[1][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) + D[1][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[2][ele[0]] += -D[2][0] * (eleU[1] - eleU[0])/len - D[2][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) - D[2][2] * (elePhi[1] - elePhi[0]) / len;
	    	ESED[2][ele[1]] += -D[2][0] * (eleU[1] - eleU[0])/len - D[2][1] * (eleU[1] + eleU[0])/ (x[0] + x[1]) - D[2][2] * (elePhi[1] - elePhi[0]) / len;
	    	
	    	count[ele[0]] += 1.0;
	    	count[ele[1]] += 1.0;
	    }
	    
	    for(int i = 0, n = mesh.getNv(); i < n; i++) {
	    	ESED[0][i] /= count[i];
	    	ESED[1][i] /= count[i];
	    	ESED[2][i] /= count[i];
	    }
	}
}
